home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Analog / Fourier.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  22.8 KB  |  746 lines
  1. (*  :Title:    Continuous Time Fourier Transform  *)
  2.  
  3. (*  :Authors:    Brian Evans, Wallace McClure, James McClellan  *)
  4.  
  5. (*  :Summary:    *)
  6.  
  7. (*  :Context:    SignalProcessing`Analog`Fourier`  *)
  8.  
  9. (*  :PackageVersion:  2.7    *)
  10.  
  11. (*
  12.     :Copyright:    Copyright 1990-1991 by Brian L. Evans
  13.         Georgia Tech Research Corporation
  14.  
  15.     Permission to use, copy, modify, and distribute this software
  16.     and its documentation for any purpose and without fee is
  17.     hereby granted, provided that the above copyright notice
  18.     appear in all copies and that both that copyright notice and
  19.     this permission notice appear in supporting documentation,
  20.     and that the name of the Georgia Tech Research Corporation,
  21.     Georgia Tech, or Georgia Institute of Technology not be used
  22.     in advertising or publicity pertaining to distribution of the
  23.     software without specific, written prior permission.  Georgia
  24.     Tech makes no representations about the suitability of this
  25.     software for any purpose.  It is provided "as is" without
  26.     express or implied warranty.
  27.  *)
  28.  
  29. (*
  30.     :History:    start -- February 22, 1990
  31.         made into package -- April 18, 1990
  32.         added Dialogue ability -- October, 1990
  33.         added many more transform pairs -- February, 1991
  34.         incorporated duality so only one set of transform
  35.           pairs, properties, etc., is encoded --  December, 1991
  36.  *)
  37.  
  38. (*  :Keywords:    Fourier transform, frequency response  *)
  39.  
  40. (*
  41.     :Source:    {Signals and Systems} by A. V. Oppenheim and A. S. Willsky.
  42.           Prentice Hall.
  43.         {The Fourier Transform and Its Applications} by R. Bracewell
  44.  *)
  45.  
  46. (*  :Warning:    *)
  47.  
  48. (*  :Mathematica Version:  1.2 or 2.0  *)
  49.  
  50. (*
  51.     :Limitation:   If the function to  be transformed does not match either
  52.            a specific Fourier transform or the LaPlace transform, a 
  53.            nasty expression is returned.
  54.  *)
  55.  
  56. (*
  57.     :Discussion: When loaded into Mathematica,  this package will allow the 
  58.          computation of  continuous  Fourier  and  inverse  Fourier
  59.          transforms.   It  turns  out  that in order to show a user
  60.          each step  of  transformation process, the transform rules
  61.          needed to be placed in a list of deferred (delayed) rules.
  62.  
  63.     MyFourier[f, t, w, s] computes  the  one-dimensional Fourier trans-
  64.     form  of  f  using  t  as  the time variable and w as the frequency
  65.     variable.
  66.  
  67.     MyInverseCTFT[f, w, t, s] computes  the   inverse   one-dimensional
  68.     continuous-time Fourier transform using w as the frequency variable
  69.     and t as the time variable.
  70.  
  71.     The Fourier  transform  formulas  used  follow Oppenheim & Willsky.
  72.     To convert Bracewell's formulas,  which use 2 Pi s x instead of w t
  73.     in the exponential function,
  74.  
  75.       (1) use  x = t  and  s = w / (2 Pi)  OR
  76.       (2) use  s = w  and  x = t / (2 Pi)  and  divide the time
  77.           function by 2 Pi.
  78.           
  79.     Note that in (2), dividing the time function by 2 Pi is equivalent
  80.     to multiplying the frequency response by 2 Pi.   Either way places
  81.     the 1/(2 Pi) normalization term in front of the inverse formula.
  82.  *)
  83.  
  84. (*  :Functions:    CTFTransform InvCTFTransform  *)
  85.  
  86.  
  87.  
  88. (*  B E G I N     P A C K A G E  *)
  89.  
  90. BeginPackage[ "SignalProcessing`Analog`Fourier`",
  91.           "SignalProcessing`Analog`LaPlace`",
  92.           "SignalProcessing`Analog`InvLaPlace`",
  93.           "SignalProcessing`Analog`LSupport`",
  94.           "SignalProcessing`Digital`DSupport`",
  95.           "SignalProcessing`Support`TransSupport`",
  96.           "SignalProcessing`Support`ROC`",
  97.           "SignalProcessing`Support`SigProc`",
  98.           "SignalProcessing`Support`SupCode`" ]
  99.  
  100.  
  101. If [ TrueQ[ $VersionNumber >= 2.0 ],
  102.      Off[ General::spell ];
  103.      Off[ General::spell1 ] ];
  104.  
  105.  
  106. (*  U S A G E     I N F O R M A T I O N  *)
  107.  
  108. CTFTData::usage =
  109.     "CTFTData is the data-tag for the CTFT."
  110.  
  111. CTFTransform::usage =
  112.     "CTFTransform[f, t] or CTFTransform[f, t, w] computes the \
  113.     continuous-time Fourier transform of f. \
  114.     It has time argument(s) t and frequency argument(s) w."
  115.  
  116. InvCTFTransform::usage =
  117.     "InvCTFTransform[f, w, t] computes the inverse continuous time \
  118.     Fourier transform of f with frequency arguments w and time \
  119.     arguments t."
  120.  
  121. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  122.  
  123.  
  124. Begin[ "`Private`" ]
  125.  
  126.  
  127. (*  Messages  *)
  128.  
  129. Transform::laplace = 
  130.     "The Region of Convergence in the Laplace transform `` does not contain the imaginary axis so it cannot be converted into a Fourier transform."
  131.  
  132.  
  133. (*  Options for the forward and inverse Fourier transforms  *)
  134.  
  135. displayFixUp[ op_ ] :=
  136.     Append[ fixUp[op], TransformLookup -> Replace[TransformLookup, op] ]
  137.  
  138. CTFTransform/: Options[ CTFTransform ] :=
  139.     displayFixUp[ Join[ { Definition -> False, Dialogue -> False },
  140.                 Options[ InvLaPlace ] ] ]
  141.  
  142. InvCTFTransform/: Options[ InvCTFTransform ] :=
  143.     displayFixUp[ Join[ { Definition -> False,
  144.                   Dialogue -> False,
  145.                   Terms -> False},
  146.                 Options[ InvLaPlace ] ] ]
  147.  
  148.  
  149. (*  Printed forms of local objects so that Dialogue -> All works neatly  *)
  150.  
  151. Format[ dualFT[ x_, w_, rest__, InvCTFTransform ] ] :=
  152.     SequenceForm[ ColumnForm[{"F" Superscript[-1],
  153.                   "  " ~StringJoin~ ToString[w]}],
  154.               { x } ]
  155.  
  156. Format[ dualFT[ x_, t_, rest__, CTFTransform ] ] :=
  157.     SequenceForm[ ColumnForm[{"F",
  158.                   "  " ~StringJoin~ ToString[t]}],
  159.               { x } ]
  160.  
  161. Format[ derivative[x_, w_Symbol, k_] ] := 
  162.     SequenceForm[ ColumnForm[{"D" Superscript[k],
  163.                   "  " ~StringJoin~ ToString[w]}],
  164.               { x } ]
  165.  
  166. Format[ integrate[x_, limits_] ] := 
  167.     SequenceForm[ "Integrate[", x, ", ", limits, "]" ]
  168.  
  169. Format[    substitute[x_, w_, wnew_] ] :=
  170.     SequenceForm[ {x}, Subscript[w -> wnew] ]
  171.  
  172. Format[    CTFTtoL[x_, w_, s_] ] :=
  173.     SequenceForm[ {x}, Subscript[w -> s / I] ]
  174.  
  175. Format[ SignalProcessing`Analog`Fourier`Private`s ] := "s"
  176.  
  177.  
  178. (*  State definition for the unified rule base *)
  179.  
  180. expandfield = 1        (* apply expand to the expression        *)
  181. forlaplacefield = 2    (* all else has failed so try forward Laplace    *)
  182. invlaplacefield = 3    (* all else has failed so try inverse Laplace    *)
  183. definitionfield = 4    (* apply the Fourier transform definition    *)
  184.  
  185. statevariables = 4
  186.  
  187. (*  Global variables  *)
  188. origVars = {}
  189. tranVars = {}
  190.  
  191. (*  Supporting routines  *)
  192.  
  193. explain[ dualFT[x_, t_, rest__, CTFTransform]] :=
  194.     Message[ Transform::incomplete, "forward Fourier transform", x, t ]
  195. explain[ dualFT[x_, w_, rest__, InvCTFTransform]] :=
  196.     Message[ Transform::incomplete, "inverse Fourier transform", x, w ]
  197.  
  198. fixUp[ op_ ] := { Apart -> Replace[Apart, op],
  199.           Definition -> Replace[Definition, op],
  200.           Dialogue -> Replace[Dialogue, op],
  201.           Simplify -> Replace[Simplify, op],
  202.           Terms -> Replace[Terms, op] }
  203.  
  204. initstate[] := Table[True, {statevariables}]
  205. initstate[ CTFTransform ] :=
  206.     SetStateField[ initstate[], invlaplacefield, False ]
  207. initstate[ InvCTFTransform ] :=
  208.     SetStateField[ initstate[], forlaplacefield, False ]
  209.  
  210. LtoCTFT[ SignalProcessing`Analog`LaPlace`Private`mylaplace[x_, lrest__], t_, w_, s_, rest__ ] :=
  211.     CTFTtoL[ dualFT[x, t, w, s, rest], w, s ]
  212. LtoCTFT[ x_, t_, w_, s_, trest__ ] := x
  213.  
  214. mustbenegative[x_] := ! TrueQ[ Positive[x] ]
  215. mustbepositive[x_] := ! TrueQ[ Negative[x] ]
  216.  
  217. posttransform[x_] := ( x /. postrules )
  218. postrules = {
  219.     derivative[x_, t_, 1] :>
  220.         D[ posttransform[x], t ],
  221.     derivative[x_, t_, k_Integer] :>
  222.         D[ posttransform[x], {t, k} ],
  223.     integrate[x_, limits_] :>
  224.         Integrate[ Expand[posttransform[x]], limits ],
  225.     substitute[x_, -v_, v_] :>
  226.         ( posttransform[x] /. v -> -v ),
  227.     substitute[x_, t_, tnew_] :>
  228.         ( posttransform[x] /. t -> tnew ),
  229.     CTFTtoL[x_] :>
  230.         Transform[ x, -Infinity, Infinity ],
  231.     CTFTtoL[x_, w_, s_] :>
  232.         Transform[ substitute[x, w, s / I], -Infinity, Infinity ],
  233.     (a_ + b_. Sign[x_]) :> 2 a CStep[-x] /; ( a == -b ),
  234.     (a_ + a_. Sign[x_]) :> 2 a CStep[x]
  235. }
  236.  
  237. ROCincludesImagAxis[ trans_, op_ ] :=
  238.     Block [    {cond, rm, rp, zerocond},
  239.         rm = GetRMinus[trans];
  240.         rp = GetRPlus[trans];
  241.         cond = ( Negative[rm] && Positive[rp] );
  242.         zerocond = ( ZeroQ[rm] || ZeroQ[rp] );
  243.         If [ zerocond && SameQ[ Replace[Dialogue, op], All ],
  244.              Print[ "assuming that the origin which is an endpoint" ];
  245.              Print[ "of the region of convergence (ROC) of the" ];
  246.              Print[ "Laplace transform is included in the ROC" ] ];
  247.         If [ cond || zerocond,
  248.              True,
  249.              False,
  250.              Assuming[cond, op]; True ] ]
  251.  
  252. transform[dualFT[x_, rest__]] := Replace[dualFT[x, rest], DualFTRules]
  253. transform[x_] := x
  254.  
  255. validvarQ[ x_Symbol ] := True
  256. validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
  257. validvarQ[ x_ ] := False
  258.  
  259.  
  260. (*  Extension of Fourier transform operators  *)
  261.  
  262. Unprotect[FT, InvFT]
  263. FT/: TheFunction[ FT[t_, w_] [f_] ] := CTFTransform[f, t, w]
  264. InvFT/: TheFunction[ InvFT[w_, t_] [f_] ] := InvCTFTransform[f, w, t]
  265. Protect[FT, InvFT]
  266.  
  267.  
  268. (*  Extension of TheFunction operator  *)
  269.  
  270. CTFTData/: TheFunction[ CTFTData[x_, var_] ] := x
  271.  
  272.  
  273. (*  Magnitude/Phase plot of a CTFT object *)
  274.  
  275. mppoptions = {    Domain -> Continuous, DomainScale -> Linear,
  276.         MagRangeScale -> Log, PhaseRangeScale -> Degree,
  277.         PlotRange -> All }
  278.  
  279. CTFTData/: MagPhasePlot[ CTFTData[trans_, vars_], rest___ ] :=
  280.     Apply [    MagPhasePlot, { trans, rest } ~Join~ mppoptions ]
  281.  
  282.  
  283. (*  One-dimensional Fourier transform driver  *)
  284.  
  285. DualFT1d[ f_, t_, w_, s_, tlist_, wlist_, oplist_, dir_ ] :=
  286.     Block [    {different, newexpr, newrules, oldexpr, simops, state, trace},
  287.  
  288.         Off[ Power::infy ];
  289.         trace = SameQ[ Replace[Dialogue, oplist], All ];
  290.         state = initstate[dir];
  291.  
  292.         (*  Generate the 1-D FT rules  *)
  293.         newrules = TransformFixUp[ t, w, oplist, dualFT, False,
  294.                        dir, Null, Null ];
  295.                 If [ SameQ[dir, InvCTFTransform],
  296.              newrules = newrules /. w -> -w ];
  297.         DualFTRules = Join[ newrules, OriginalDualFTRules ];
  298.  
  299.         (*  Embed the function in the dual ctft object  *)
  300.         If [ trace && SameQ[dir, InvCTFTransform],
  301.              Print[ "The Duality Theorem will be used to" ];
  302.              Print[ "find the inverse Fourier transform:" ] ];
  303.         oldexpr = dualFT[f, t, w, s, fixUp[oplist], state, dir];
  304.         If [ SameQ[dir, InvCTFTransform],
  305.              oldexpr = dualFT[ f, t, w, s, fixUp[oplist],
  306.                             state, CTFTransform ];
  307.              oldexpr = substitute[ oldexpr, w, -w ] / (2 Pi);
  308.              If [ trace, Print[oldexpr]; Print["which becomes"] ] ];
  309.  
  310.         (*  Simplify the function before transformation  *)
  311.  
  312.         simops = { Dialogue -> trace, Variables -> tlist };
  313.         If [ TrueQ[ $VersionNumber >= 2.0 ],
  314.              PrependTo[simops, Trig -> False] ];
  315.         newexpr = SPSimplify[ oldexpr, simops ];
  316.  
  317.         (*  Transform the function  *)
  318.         different = True;
  319.         If [ trace, Print[newexpr]; Print["which becomes"] ];
  320.         While [ different,
  321.             oldexpr = newexpr;
  322.             newexpr = oldexpr /. ( dualFT[a__] :>
  323.                   Replace[dualFT[a], DualFTRules] );
  324.             different = ! SameQ[newexpr, oldexpr];
  325.             If [ different && trace,
  326.                  Print[newexpr]; Print["which becomes"] ] ];
  327.  
  328.         (*  Fix up the transform  *)
  329.         different = True;
  330.         While [ different,
  331.             oldexpr = newexpr;
  332.             newexpr = posttransform[oldexpr];
  333.             different = ! SameQ[newexpr, oldexpr];
  334.             If [ different && trace,
  335.                  Print[newexpr]; Print["which becomes"] ] ];
  336.  
  337.         (*  Check for invalid Laplace transform substitutions  *)
  338.         newexpr /. LTransData[x_, rest__] :>
  339.                      Message[ Transform::laplace, x ];
  340.  
  341.         (*  Simplify the resulting transform  *)
  342.         If [ trace, Print[newexpr] ];
  343.         newexpr = SPSimplify[ newexpr,
  344.                       Variables -> wlist,
  345.                       Dialogue -> trace ];
  346.  
  347.         On[ Power::infy ];
  348.  
  349.         newexpr ]
  350.  
  351.  
  352.  
  353.  
  354. (*  B E G I N     F O R W A R D     F O U R I E R     T R A N S F O R M  *)
  355.  
  356.  
  357. (*  CTFTransform  *)
  358. CTFTransform[ f_ ] :=
  359.     Message[ Transform::novariables, "t (time)", GetVariables[f] ]
  360. CTFTransform[ f_, t_ ] :=
  361.     CTFTransform[ f, t, DummyVariables[Length[t], Global`w] ]
  362. CTFTransform[ f_, t_, w_, options___ ] :=
  363.     Message[ Transform::badvar, "time", CTFTransform, t ] /;
  364.     ! validvarQ[t]
  365. CTFTransform[ f_, t_, w_, options___ ] :=
  366.     Message[ Transform::badvar, "frequency", CTFTransform, w ] /;
  367.     ! validvarQ[w]
  368. CTFTransform[ f_, t_, w_, options___ ] :=
  369.     Block[ {op, svars, trans, vars},
  370.         vars = Which [ Length[t] > Length[w],
  371.                  DummyVariables[Length[t], Global`w],
  372.                    Length[t] < Length[w],
  373.                  Take[w, Length[t]],
  374.                    True,
  375.                  w ];
  376.         svars = DummyVariables[Length[t], s];
  377.  
  378.         origVars = t;
  379.         tranVars = vars;
  380.         op = ToList[options] ~Join~ Options[CTFTransform];
  381.         trans = If [ Length[t] > 1, 
  382.                  MultiDCTFT[f, t, vars, svars, op],
  383.                  MyFourier[f, t, vars, svars, op] ];     
  384.  
  385.         If [ InformUserQ[op],
  386.              Scan[explain, trans, Infinity] ];
  387.  
  388.         CTFTData[ FourierSimplify[trans], FVariables[vars] ] ]
  389.  
  390. (*  Multidimensional continuous-time Fourier transform  *)
  391. MultiDCTFT[f1_, t_, w_, s_, op_] :=
  392.     Block [    {i, func = f1, length},
  393.         length = Length[t];
  394.         For [ i = 1, i <= length, i++,
  395.               func = MyFourier[ func, t[[i]], w[[i]],
  396.                     s[[i]], t, w, op ] ];
  397.         func ]
  398.  
  399. (*  One-dimensional continuous-time Fourier transform  *)
  400. MyFourier[ f_, t_, w_, s_, oplist_ ] :=
  401.     DualFT1d[ f, t, w, s, {t}, {w}, oplist, CTFTransform ]
  402.  
  403. MyFourier[ f_, t_, w_, s_, tlist_, wlist_, oplist_ ] :=
  404.     DualFT1d[ f, t, w, s, tlist, wlist, oplist, CTFTransform ]
  405.  
  406.  
  407. (*  E N D     F O R W A R D     F O U R I E R     T R A N S F O R M  *)
  408.  
  409.  
  410.  
  411.  
  412. (*  B E G I N     I N V E R S E     F O U R I E R     T R A N S F O R M  *)
  413.  
  414.  
  415. (*  InvCTFTransform  *)
  416. InvCTFTransform[ CTFTData[x_, FVariables[w_]] ] := InvCTFTransform[x, w]
  417.  
  418. InvCTFTransform[ CTFTData[x_, FVariables[w_]], rest__ ] :=
  419.     InvCTFTransform[x, rest]
  420.  
  421. InvCTFTransform[ x_, w_ ] :=
  422.     InvCTFTransform[ x, w, DummyVariables[Length[w], Global`t] ]
  423.  
  424. InvCTFTransform[ f_, w_, t_, options___ ] :=
  425.     Message[ Transform::badvar, "frequency", InvCTFTransform, w ] /;
  426.     ! validvarQ[w]
  427.  
  428. InvCTFTransform[ f_, w_, t_, options___ ] :=
  429.     Message[ Transform::badvar, "time", InvCTFTransform, t ] /;
  430.     ! validvarQ[t]
  431.  
  432. InvCTFTransform[ f_, w_, t_, options___ ] :=
  433.     Block[ {op, s, svars, trans, vars},
  434.         op = ToList[options] ~Join~ Options[InvCTFTransform];
  435.  
  436.         vars = Which [ Length[w] > Length[t],
  437.                  DummyVariables[Length[w], Global`t],
  438.                    Length[w] < Length[t],
  439.                  Take[t, Length[w]],
  440.                    True,
  441.                  t ];
  442.  
  443.         svars = DummyVariables[Length[w], s];
  444.         origVars = w;
  445.         tranVars = vars;
  446.         trans = If [ Length[w] > 1, 
  447.                  MultiDInvCTFT[f, w, vars, svars, op],
  448.                  MyInverseCTFT[f, w, vars, svars, op] ];
  449.  
  450.         If [ InformUserQ[op],
  451.              Scan[explain, trans, Infinity] ];
  452.  
  453.         FourierSimplify[trans] ]
  454.  
  455. (*  MultiDInvCTFT  *)
  456. MultiDInvCTFT[ f1_, w_, t_, s_, op_ ]:= 
  457.     Block [    {i, func = f1, length},
  458.         length = Length[t];
  459.         For [ i = 1, i <= length, i++,
  460.               func = MyInverseCTFT[ func, w[[i]], t[[i]],
  461.                         s[[i]], w, t, op ] ];
  462.         func ]
  463.  
  464. (*  driver for 1-D inverse CTFT rules  *)
  465. MyInverseCTFT[ f_, w_, t_, s_, op_ ] :=
  466.     DualFT1d[ f, w, t, s, {w}, {t}, op, InvCTFTransform ]
  467.  
  468. MyInverseCTFT[ f_, w_, t_, s_, wlist_, tlist_, op_ ] :=
  469.     DualFT1d[ f, w, t, s, wlist, tlist, op, InvCTFTransform ]
  470.  
  471.  
  472.  
  473.  
  474.  
  475. (*  E N D     I N V E R S E     F O U R I E R     T R A N S F O R M  *)
  476.  
  477.  
  478.  
  479.  
  480. (*  B E G I N     D U A L     R U L E     B A S E  *)
  481.  
  482. DualFTRules = {}
  483.  
  484. OriginalDualFTRules = {
  485.  
  486. (*      Clean up any calls to the Laplace transform rule bases          *)
  487. (*      Must go before any other rules to avoid spurious Delta functions  *)
  488.  
  489. dualFT[ LTransData[x_, rest__], t_, w_, s_, op_, __ ] :>
  490.     If [ ROCincludesImagAxis[ LTransData[x, rest], op ],
  491.          x /. s -> I w,
  492.          LTransData[x, rest] ] ,
  493.  
  494. dualFT[ SignalProcessing`Analog`InvLaPlace`Private`myinvlaplace[x_,
  495.         s_, restL__], t_, w_, s_, rest__ ] :>
  496.     dualFT[ x /. s -> I t, t, w, s, rest ],
  497.  
  498.  
  499. (*  A.  Transform pairs            *)
  500. (*     1.  constant functions        *)
  501. dualFT[ 0, t_, w_, __ ] :> 0,
  502. dualFT[ c_, t_, w_, __ ] :> 2 Pi c Delta[w] /; FreeQ[c, t],
  503.  
  504. (*     2.  two delta functions    *)
  505. (*         c delta(t + a) + d delta(t + b) where c = -d and a = -b    *)
  506. (*         c delta(t - b) - c delta(t + b) <---> -2 I c sin(b w)     *)
  507. dualFT[ c_. Delta[t_ + a_] + d_. Delta[t_ + b_], t_, w_, __ ] :>
  508.     -2 I c Sin[b w] /;
  509.     ( a == -b ) && ( c == 1 || c == -1 ) && ( c == -d ) && FreeQ[{a, c}, t],
  510.  
  511. (*         c delta(t + a) + c delta(t + b) where c = -d and a = -b    *)
  512. (*         c delta(t - b) + c delta(t + b) <---> 2 c cos(b w)     *)
  513. dualFT[ c_. Delta[t_ + a_] + c_. Delta[t_ + b_], t_, w_, __ ] :>
  514.     2 c Cos[b w] /;
  515.     ( a == -b ) && FreeQ[{a, c}, t],
  516.  
  517. (*     3.  one delta function    *)
  518. dualFT[ c_. Delta[t_ + t0_.], t_, w_, __ ] :>
  519.     c Exp[I t0 w] /;
  520.     FreeQ[{c, t0}, t],
  521.  
  522. (*     4.  impulse train        *)
  523. dualFT[ x_. Periodic[T_, t_][y_. Delta[t_ + a_.]], t_, w_, rest__ ] :> 
  524.     (1 / Abs[T]) Exp[I a w] *
  525.       Periodic[2 Pi / T, w][ dualFT[x y, t, w, rest] ] /;
  526.     FreeQ[{a, n, y, T}, t],
  527.  
  528. dualFT[ x_. Summation[n_Symbol, -Infinity, Infinity, 1][
  529.               y_. Delta[t_ + a_. + n_ T_.]],
  530.     t_, w_, rest__ ] :> 
  531.     (1 / Abs[T]) Exp[I a w] *
  532.       Summation[n, -Infinity, Infinity, 1][
  533.         substitute[ dualFT[x y, t, w, rest], w, w + 2 Pi n / T ] ] /;
  534.     FreeQ[{a, n, T}, t] && FreeQ[y, n],
  535.  
  536. (*     5.  pulse functions        *)
  537. dualFT[ CStep[t_ + a_] - CStep[t_ + b_] + x_., t_, w_, rest__ ] :>
  538.     2 a Sinc[a w] + dualFT[x, t, w, rest] /;
  539.     FreeQ[a, t] && (a == -b),
  540.  
  541. dualFT[ x_. + c_. ( CStep[t_ + a_.] - CStep[t_ + b_.] ), t_, w_, rest__ ] :>
  542.     dualFT[c CPulse[Abs[b - a], t - Min[-a, -b]], t, w, rest] +
  543.       dualFT[x, t, w, rest] /;
  544.     FreeQ[{a, b}, t],
  545.  
  546. dualFT[ CPulse[T_, t_ + a_.], t_, w_, __ ] :>
  547.     T Exp[-I w T / 2] Exp[I a w] Sinc[w T / 2] /;
  548.     FreeQ[{a, T}, t],
  549.  
  550. (*     6.  two-sided Gaussian        *)
  551. dualFT[ Exp[a_. t_^2], t_, w_, s_, op_, rest___ ] :>
  552.     Block [ {},
  553.         Assuming[ Negative[a], op ];
  554.         Sqrt[Pi] Exp[w^2 / (4 a)] / Sqrt[-a] ] /;
  555.     FreeQ[a, t] && mustbenegative[a],
  556.  
  557. (*     7.  cosine x over x in time    *)
  558. dualFT[ Cos[a_. t_] / t_, t_, w_, __ ] :>
  559.     Abs[a] I ( CStep[-Abs[a] - w] - CStep[w - Abs[a]] ) /;
  560.     FreeQ[a, t],
  561.  
  562. (*     8.  two-sided cosine        *)
  563. dualFT[ Cos[w0_. t_ + b_.], t_, w_, __ ] :>
  564.     Exp[I w b / w0] Pi ( Delta[w - w0] + Delta[w + w0] ) /;
  565.     FreeQ[{b, w0}, t],
  566.  
  567. (*     9.  two-sided sine        *)
  568. dualFT[ Sin[w0_. t_ + b_.], t_, w_, __ ] :>
  569.     Exp[I w b / w0] Pi / I ( Delta[w - w0] - Delta[w + w0] ) /;
  570.     FreeQ[{b, w0}, t],
  571.  
  572. (*    10.  two-sided sinc        *)
  573. dualFT[ Sin[a_. t_] / t_, t_, w_, __ ] :> 
  574.     Pi CPulse[2 a, w + a] /;
  575.     FreeQ[a, t],
  576.  
  577. dualFT[ Sinc[a_. t_ + b_.], t_, w_, __ ] :> 
  578.     Exp[I w b / a] Pi CPulse[2 a, w + a] / a /;
  579.     FreeQ[{a,b}, t],    (* CPulse[2a, w+a] = CStep[w+a] - CStep[w-a] *)
  580.  
  581. dualFT[ (Sinc[a_. t_ + b_.])^2, t_, w_, __ ] :> 
  582.     Exp[I w b / a] (Pi / a)^2 *
  583.       ( Unit[-2][w + 2 a] - 2 Unit[-2][w] + Unit[-2][w - 2 a] ) /;
  584.     FreeQ[{a,b}, t],
  585.  
  586. (*    11.  other two-sided forms    *)
  587. dualFT[ Exp[-Abs[a_. t_]] Sinc[a_. t_], t_, w_, __ ] :>
  588.     ArcTan[ 2 a^2 / w^2 ] / a /;
  589.     FreeQ[a, t],
  590.  
  591. dualFT[ Sign[t_ + t0_.], t_, w_, __ ] :>
  592.     2 I Exp[I t0 w] / w /;
  593.     FreeQ[t0, w],
  594.  
  595. dualFT[ 1 / t_, t_, w_, __ ] :> I Pi Sign[w],
  596.  
  597. dualFT[ BesselJ[0, a_. t_], t_, w_, __ ] :>
  598.     CPulse[2 a, w + a] / ( Pi Sqrt[a^2 - w^2] ) /;
  599.     FreeQ[a, t],
  600.  
  601. dualFT[ CStep[t_ + t0_.], t_, w_, __ ] :>
  602.     I Exp[I t0 w] / w  +  Pi Delta[w] /;
  603.     FreeQ[t0, t],    (* because Pi Exp[I t0 w] Delta[w] <--> Pi Delta[w] *)
  604.  
  605. dualFT[ CStep[-t_ + t0_.], t_, w_, __ ] :>
  606.     -I Exp[-I t0 w] / w  +  Pi Delta[w] /;
  607.     FreeQ[t0, t],    
  608.  
  609. (*  B.  Properties  *)
  610. (*    1.  homogeneity        *)
  611. dualFT[ a_ x_, t_, w_, rest__ ] :> 
  612.     a dualFT[ x, t, w, rest ] /; FreeQ[a, t],
  613.  
  614. (*    2.  additivity        *)
  615. dualFT[ a_ + b_, t_, w_, rest__ ] :> 
  616.     dualFT[ a, t, w, rest ] + dualFT[ b, t, w, rest ],
  617.  
  618. dualFT[ Summation[i_, il_, iu_, inc_][f_], t_, w_, rest___ ] :>
  619.     Summation[i, il, iu, inc][ dualFT[f, t, w, rest] ] /;
  620.     FreeQ[{i, il, iu, inc}, t],
  621.  
  622. (*    3.  one-sided transform    *)
  623. dualFT[ x_. CStep[t_ + t0_], t_, w_, rest__ ] :>
  624.     Exp[- I w t0] dualFT[ ( x /. t -> (t - t0)) CStep[t], t, w, rest ] /;
  625.     FreeQ[t0, t],
  626.  
  627. (*    4.  derivative        *)
  628. dualFT[ Derivative[k_][x_][t_], t_, w_, rest__ ] :>
  629.     (I w)^k dualFT[ x[t], t, w, rest ] /;
  630.     FreeQ[k, t],
  631.  
  632. (*    5.  integration        *)
  633. dualFT[ Integrate[x_, {tau_, -Infinity, t_}], t_, w_, rest__ ] :>
  634.     Block [    {trans},
  635.         trans = dualFT[x, tau, w, rest];
  636.         trans ( 1/(I w) + Pi Delta[w] ) ],
  637.  
  638. (*    6.  multiply by time    *)
  639. dualFT[ t_^a_. x_, t_, w_, s_, op_, rest___ ] :> 
  640.     Block [    {dialogue},
  641.         dialogue = SameQ[ Replace[Dialogue, op], All ];
  642.         Assuming[ Positive[a], dialogue ];
  643.         If [ dialogue, Print[ "where ", a, " is an integer" ] ];
  644.         I^a derivative[ dualFT[x, t, w, s, op, rest], w, a ] ] /;
  645.     FreeQ[a, t] && Implies[ NumberQ[a], IntegerQ[a] && a > 0 ],
  646.  
  647. (*    7.  divide by time      *)
  648. dualFT[ x_ / t_^k_., t_, w_, rest__ ] :> 
  649.     -I integrate[ dualFT[x / t^(k-1), t, w, rest], {w, -Infinity, w} ] +
  650.       I Pi Limit[x, t -> 0] /;
  651.     N[ k >= 1 ] && ! InfinityQ[ Limit[x, t -> 0] ],
  652.  
  653. (*    8.  exponential raised to an imaginary power times a function  *)
  654. dualFT[ x_. Exp[ Complex[0, w0_] w1_. t_ + restexp_. ], t_, w_, rest__ ] :>
  655.     substitute[ dualFT[x Exp[restexp], t, w, rest], w, w - w0 w1 ] /;
  656.     FreeQ[{w0, w1}, t],
  657.  
  658. (*     9.  convolution        *)
  659. dualFT[ CConvolve[t_][x_, restconv__], t_, w_, rest__ ] :> 
  660.     Apply[ Times, Map[ dualFT[#1, t, w, rest]&, { x, restconv } ] ],
  661.  
  662.  
  663. (*  C.    Strategies        *)
  664. (*    1.  expand expression    *)
  665. dualFT[ x_, t_, w_, s_, oplist_, state_, rest___ ] :>
  666.     dualFT [ Expand[x], t, w, s, oplist,
  667.          SetStateField[state, expandfield, False], rest ] /;
  668.     GetStateField[state, expandfield],
  669.  
  670. (*    2.  function times a window    *)
  671. dualFT[ x_ CPulse[l_, t_ + a_.], t_, w_, __ ] :>
  672.     Integrate[ x Exp[- I w t], {t, -a, -a + l} ] /;
  673.     FreeQ[{a, l}, t],
  674.  
  675. (*    3.  Apply the definition if Definition options is enabled  *)
  676. dualFT[ x_, t_, w_, s_, op_, state_, dir_ ] :>
  677.     Block [    {newstate, trans},
  678.         newstate = SetStateField[state, definitionfield, False];
  679.         trans = Integrate[x Exp[- I w t],
  680.                   {t, -Infinity, Infinity}];
  681.         If [ SameQ[Head[trans], Integrate],
  682.              dualFT[x, t, w, s, op, newstate, dir],
  683.              TransformDialogue[ Definition, op, Integrate,
  684.                     x, trans ] ] ] /;
  685.     GetStateField[state, definitionfield] && Replace[Definition, op],
  686.  
  687. (*    4.  Try a forward or inverse Laplace transform if all else fails  *)
  688. dualFT[ x_, t_, w_, s_, oplist_, state_, dir_ ] :> 
  689.     Block [    {ltrans, newstate},
  690.         ltrans = LaPlace[ x, t, s,
  691.                   { Dialogue -> False } ~Join~ oplist ];
  692.         If [ InformUserQ[oplist],
  693.              Print[ "( after taking the Laplace transform of " ];
  694.              Print[ "  ", x ];
  695.              Print[ "  with respect to ", t, " and" ];
  696.              Print[ "  adjusting the result )" ] ];
  697.         newstate = SetStateField[state, forlaplacefield, False];
  698.         If [ SameQ[Head[ltrans], LTransData],
  699.              If [ ROCincludesImagAxis[ ltrans, op ],
  700.               ( TheFunction[ltrans] /. s -> I w ),
  701.                   dualFT[ x, t, w, s, oplist, newstate, dir ] ],
  702.              MapAll[ LtoCTFT[#1, t, w, s, oplist, newstate, dir]&,
  703.                  ltrans ] ] ] /;
  704.     GetStateField[state, forlaplacefield],
  705.  
  706. dualFT[ x_, w_, t_, s_, op_, state_, dir_ ] :>
  707.     Block [    {invltrans, newstate, newx},
  708.         newx = (x /. w -> (s / I));
  709.         invltrans = InvLaPlace[ 2 Pi newx, s, t,
  710.                     { Dialogue -> False } ~Join~ op ];
  711.         If [ InformUserQ[op],
  712.              Print[ "( after taking the inverse Laplace transform of" ];
  713.              Print[ "  ", newx /. s -> "s" ];
  714.              Print[ "  which yields" ];
  715.              Print[ "  ", invltrans ];
  716.              Print[ "  and adjusting the result )" ] ];
  717.         newstate = SetStateField[state, invlaplacefield, False];
  718.         If [ FreeQ[invltrans, s],
  719.              substitute[invltrans, t, -t],   (* FT uses L^-1 so t=-t *)
  720.              dualFT[ x, w, t, s, op, newstate, dir ] ] ] /;
  721.     GetStateField[state, invlaplacefield]
  722.  
  723. }
  724.  
  725.  
  726. (*  E N D     P A C K A G E  *)
  727.  
  728. End[]
  729. EndPackage[]
  730.  
  731. If [ TrueQ[ $VersionNumber >= 2.0 ],
  732.      On[ General::spell ];
  733.      On[ General::spell1 ] ];
  734.  
  735.  
  736. (*  H E L P     I N F O R M A T I O N  *)
  737.  
  738. Combine[ SPfunctions, { CTFTransform, InvCTFTransform } ]
  739. Protect[ CTFTransform, InvCTFTransform ]
  740.  
  741.  
  742. (*  E N D I N G     M E S S A G E  *)
  743.  
  744. Print[ "The continuous Fourier transform (CTFT) functions are loaded." ]
  745. Null
  746.